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The development and use of Monte Carlo algorithms plays a visible role in the study of non- 
Markovian quantum dynamics due to the provided insight and powerful numerical methods for 
solving the system dynamics. In the Markovian case, the connections between the various types 
of methods are fairly well-understood while for non-Markovian case there has so far been only a 
few studies. We focus here on two jumplike unravelings of non-Markovian dynamics, the non- 
Markovian quantum jump (NMQJ) method and the property state method by Gambetta, Askerud, 
and Wiseman (GAW). The results for simple quantum optical systems illustrate the connections 
between the realizations of the two methods and also highlight how the probability currents between 
the system and environment, or between the property states of the total system, associate to the 
decay rates of time- local master equations, and consequently to the jump rates of the NMQJ method. 

PACS numbers: 03.65.Yz, 42.50.Lc 



I. INTRODUCTION 

The theory of open quantum systems deals with the 
dynamics of the reduced system which is coupled to its 
environment pQ . This leads often to decoherence and the 
loss of quantum properties (2HS] , though there also exists 
schemes to exploit system-reservoir interaction for quan- 
tum engineering [SHE]- Recently, non-Markovian dynam- 
ics, where memory effects play a crucial role, has become 
under very active research (9lf2T]. On one hand this is 
due to the fact that the fundamental understanding of 
non-Markovianity is still missing and on the other hand 
non-Markovianity may be useful for various quantum in- 
formation or quantum engineering tasks |20|, 122] . 

The solving of non-Markovian dynamics is often a chal- 
lenging task and there exists a large number of both ana- 
lytical methods [1, 23 -26] and numerical Monte Carlo al- 
gorithms for this purpose [TTHT31 |2"TH3"T] . Roughly speak- 
ing, the Monte Carlo methods can be divided into discon- 
tinuous jumplike unravelings or continuous diffusion type 
unravelings. For the Markovian case without memory ef- 
fects, the connections between the methods are fairly well 
understood [H [381443] while the same can not be said of 
the non-Markovian methods despite of a few early studies 
[551111]. 

We focus here on two jumplike unravelings and illus- 
trate their connections by studying simple quantum op- 
tical systems. The method by Gambetta, Askerud, and 
Wiseman (GAW) is based on the generating stochastic 
realizations for the total system state vectors and mon- 
itoring the random jumps between the property states 
of the total system [35] . On the other hand, the re- 
cently developed non-Markovian quantum jump (NMQJ) 
method generates jumplikc realizations for the state vec- 
tors within the Hilbert space of the reduced system. 
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We show here that there is an inherent connection 
between the reduced system part of the GAW realiza- 
tions and the NMQJ realizations. Moreover, we also 
study how the probability currents between the property 
states are associated to the decay rates of the time-local 
master equations, and how the jump rates between the 
GAW and NMQJ methods are connected. We stress that 
the NMQJ method can be currently used for systems 
for which time local non-Markovian master equation can 
be derived whereas GAW method has greater generality. 
Our results provide new insight for non-Markovian dy- 
namics in terms of the information flow between the sys- 
tem and the reservoir, and hopefully stimulates further 
studies of connections between Monte Carlo methods for 
non-Markovian dynamics. 

The paper is organized in the following way. In Sec- 
tion [n] we describe the basic ingredients of the GAW 
method and in Sec. |III| of the NMQJ method. By study- 
ing simple quantum optical systems, in Sec. |IV| we show 
how the methods are connected and finally, Sec. [V] con- 
cludes the paper. 



II. UNRAVELING IN THE TOTAL SYSTEM 
SPACE: GAW METHOD 

The method by Gambetta, Askerud and Wiseman 
(GAW) is based on generating piecewise deterministic 
realizations, or jumplike unraveling within the Hilbert 
space of the total system, describing the discontinuous 
transitions between the property states of the system. 
We give here the basic ingredients of the method suitable 
for undriven quantum optical systems with spectral mode 
unraveling. We note that GAW method can also be ap- 
plied to driven systems and with temporal mode unravel- 
ing. More details can be found from the Refs. [3BJ H5J 06] . 

We focus on the dynamics of simple undriven quantum 
optical systems, e.g., two-level and V-systems which are 
coupled to a continuum of electromagnetic field modes 
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at zero temperature. The dynamics of the state vector 
of the system and the environment in T-Lg ® 7~ts j where 
Hg and T-Lg are the Hilbert spaces of the system and 
the environment respectively, is given by the Schrbdinger 
equation 



dt 



mt))=-iH\*(t)). 



(i) 



Here we have set h = 1. The Hamiltonian H = Hg + 
Hg + Hg£ includes the free evolution of the system Hg 
and the environment Hg, and the system-environment 
interaction Hgg. The free evolution of the n-level sys- 
tem is governed by Hg = Ylk=i w fc|^)(^l> where |fc) are 
the energy eigenstates of the system and ujk the corre- 
sponding energies. The free evolution of the TV-mode 
environment is given by He — Y^j=i v j a \ a j^ where the 
operators aj(cij) are the annihilation (creation) operators 
for jth mode of the environment. For simplicity we fo- 
cus on system-environment interactions that under the 
rotating wave approximation (RWA) include only transi- 
tion from the excited states to the unique ground state 
without any cascade structure. The general form of such 
interaction is 



N 



Hg £ = l J2J2(9k\l)(k\a] - gi\k)(l\ aj ). 



(2) 



fc>i j=i 



From now on we will work in the interaction picture H — > 
#/(£) = e i ^ Hs+H ^ t He- l{Hs+H£)t 1 where the dynamics is 
given by Eq. |lj with the Hamiltonian 

n N 

ff ; (()=^^fc|l)(fc|«je- n ^- S ^)(l| ¥ '^"), 

k>l j=l 

(3) 

where flj ^ = Vj — ujk and is the energy difference of 
the ground state (labeled with index 1) and fcth excited 
state of the system. Here g k is a frequency-dependent 
coupling constant. 

The total system state vector ^(i)) evolves according 
to the Eq. ([I]) with the Hamiltonian Let us define a 
projective operator valued measure (POVM) as 



= Ig (g) |tojv)(tojv|, 



(4) 



where m^r is a shorthand notation for arbitrary photon 
number configuration of the N environmental modes. In 
the systems we study we can have at maximum one ex- 
citation in the environment. However, GAW method in 
general is not limited only to one excitation |36j . We can 
now define property states of the total system which are 
conditioned on some particular photon number configu- 
ration of the environment. These are 



\* mN )=n mN \*V))/\fK 



-\<f>m N {t))s ® \m N ) 



(5) 



where N mN is a normalization factor. We denote the un- 
normalized property states with \^> mN (t)). From now on 
we drop the subscript TV for notational convenience and 
the simple index m refers to a particular configuration of 
iV environmental modes. 

The GAW method is a piecewise deterministic process 
(PDP) where the jumps take place between the different 
property states \^ m (t)) [SHUT). Let us define P(m,t) 
as the probability for the total system to be in the state 
\^m(t)) at time t and write the following equation of 
motion for P(m,t): 



dt 



(6) 



where J m ,k(t) is the probability current from \^k(t)} to 
|* m (i)), when J m ,k(t) > 0, and when J m ,k(t) < it 
is the probability current from from |^ m (f)) to \^>k(t)), 
i.e. to the opposite direction. We can define this in the 
following way 

J m ,k{t) = T mik {t)P{k,t)-T kim {t)P{m,t), (7) 

where T m ^(t) is the transition rate from \^k(t)) to 
|\I/ m (i)). From this definition it is clear that J m ,k(t) — 
— Jk,m(t)- Given J m ,k{t) and P(m, t) there are many pos- 
sible transition rates satisfying Eq. ([6]). One possibility 
is to use the following one [351 HZ] : 
When J mik (t) > 0, 



P(k,t) ' 
0, 



(8) 



and when J m ,k(t) < 0, 



T m: k(t) — 0, 

Tk,m(t) = — 



Jm,k{t) 

P(m,t)' 



(9) 



Since we know the total wave function of the system and 
the environment ^(f)), the probability of a given prop- 
erty state |* m (*))is P(m,t) = {9(t)\n m \^(t)). Using 
P(m, t) and Eq. with the Hamiltonian in Eq. ^ we 
obtain 



J ro , fe (t) = 2hn{(V(t)\n m H I (t)n k \*(t))} 



(10) 



The method for generating the realizations of the pro- 
cess begins with solving the total system Schrodinger 
equation followed by the calculation of the quantities 
Jm,k(t), T m> k(t) and P(k,t). For example, the proba- 
bility to have a jump between the property states of the 
total system |*k(*)) -> |* m (t)), when J m ,k(t) > 0, be- 
tween [t, t + St] is 5tT m! k{t). Then using random numbers 
we can decide whether a jump takes place or not. From 
Eqs. d8| and (10) we see that the term StT m k (t) includes 



the occupation probability of the source state k in the 
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ensemble and the rate term from k to m and together 
they give the transition rate for a single trajectory. 

After generation of the realizations, the state of the 
reduced system is 

Ps (t) =tr £ w m (t) | * m (t)) (tf m (t) | } 

m 

= ^2w m (t)\<f> m (t)){<l> m (j:)\, (11) 



where w m (t) 



M 



are the approximations for proba- 



bilities P(m,t), M is the size of the statistical ensemble 
and #(m) is the number of ensemble members in state 



III. UNRAVELING IN THE REDUCED 
SYSTEM SPACE: NMQJ METHOD 

The non-Markovian quantum jump (NMQJ) method 
is constructed as a piecewise deterministic process in the 
Hilbert space of the system Hs [TTHT3] and the key in- 
gredient is the association of negative decay rates to re- 
verse quantum jumps. The starting point is the time- 
local non-Markovian master equation, which can be de- 
rived, e.g. with time-convolutionless projection operator 
method (TCL) 1 . General form of such master equation 
(given here in the interaction picture) is 



dt^ (t) 



■■-i{H hS {t),ps(t)} 

+ ^A j (t)(c jPs {t)Cft j - \{ps{t),C)Cj}\ 
j 

(12) 

where H LS (t) = \ V . SAhC[C, is the Lamb shift 
Hamiltonian, Sj(t) is the Lamb shift rate, Aj(t) is the 
decay rate which can take negative values and operator 
Cj is the Lindblad or jump operator to channel j. The 
density matrix of the system at any point of time is de- 
composed as 

M ef f 

Ps (t) = J2P(\Mt)),t)Mt))(Mt)l (13) 

i=l 

where M is the size of the statistical ensemble of the 
unraveling, M e s is the dimension of the set of different 
states needed in the simulation (so called effective ensem- 
ble size), and P(\ipi(t)),t) is the probability of finding 
state \tpi(t)} (i(ji(t)\ in ps(t). The states in the ensemble 
evolve according to [TJ [42] 



dt 



\iPi{t)} = -iH eS (t)\iPi{t)) 



(14) 



The rate of jumps during positive decay in channel j 
from state \tpk(t)) to state \ipi(t)) with jump operator Cj 
is 



R{ k (t) = A j (t){Mt)\C}C j \i, k (t)). 
The corresponding quantum jump is given by 



\Mt)) -► \Mt)) 



(Mt)\c]c 3 \Mt)) 



(15) 



(16) 



The action of operator Cj thus means that the state 
\tpk(t)) is destroyed and the state \ipi(t)) is created in 
the statistical ensemble. 

During a negative decay probability period the jumps 
occur in reverse direction in the following sense: 



\Mt)} <- \Mt)) 



(17) 



yJmt)\C]Cj\Mt)) 
The rate of these reverse jumps is obtained from 



IV. CONNECTION BETWEEN THE GAW AND 
NMQJ UNRAVELINGS 

To make a connection between the two unravelings, 
we are interested (i) whether the reduced system part of 
the total system property state realizations of the GAW 
method have similarities with the NMQJ realizations, 
and (ii) if the jumps within the two methods occur with 
the same rates. As we will show below, the answer for 
both of these questions is positive. 

Comparing the rates, Eqs. ^ and ( [18] ), we note that 
the jump rates for the reverse probability flow and neg- 
ative decay rates, J m ^ < and Aj(t) < respectively, 
have similar structure. They both are inversely propor- 
tional to the probability to be in the source state of the 
jump. In the GAW method, the given property state is 
associated to specific mode to have the excitation (unless 
the environment is in the vacuum state). In the NMQJ 
realizations, we know whether the system or the environ- 
ment has the excitation while in the latter case we do not 
know which specific mode has the excitation. 

In order to reveal the detailed connection between the 
GAW and NMQJ methods, let us define the following 
operators 



n =I S <8 |0i, 2 , OjvXOi.Oa, 0*1 = Is ® |0> <0|, 

N N 
III =^/ S <g>4l°>v D l a fc = S J S®l 1 fc)<lfc|. 



(19) 



k=l 



k=l 



From Eq. Q we see that IIo = 7r m= o...Q and LTi = ^2 k 7Tfc, 
where fc = • • • 1& • • • (ie. k labels all one mode config- 
urations of the environment). We can now ask what is 
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the probability P(Q, t) to find zero photons at time t in 
the environment? This is given by 



p(o,i) = |n 1 *(«)>• 



(20) 



Similarly, the total probability P(l, t) of having one pho- 
ton in the environment, but not knowing in which mode, 
is 



p(i,t) = <*(t)|n 1 |*(t)). 



(21) 



The connection between the GAW and the NMQJ meth- 
ods is found by re-formulating the GAW method for the 
following combined property states: 



l*o(*)> 



1 



=n„|\Kt)>, 



=-^n 1 |*(t)>. 



(22) 



For this purpose, we must calculate the combined proba- 
bility current from the system to the environment. This 
is obtained by considering the total probability current 
from the TV-mode vacuum states to all lfc-states 



N 



•7i,o(*) = X) J w>(*)- 



(23) 



k=l 



It can be easily shown that the combined probability cur- 
rent satisfies Ji,o(t) = —Jo,i(t) and 



dt 
d_ 

d7 



P(l,t) = Ji,o{t), 
P(0,f) = -J 1 , (f). 



(24) 



We have ELl^.*) = J« P (M) and the r - h - s - of 
both equations follow from the definition of Eqs. ([6| 
and ( |23| ). The transition rates have similar structure 
as m Eqs. |[8j| and @ but we must replace probability 
with combined probability and probability current with 
combined probability current. 

As we will show below for specific examples, the GAW 
transition rates defined with combined quantities corre- 
spond to the transition rates of the NMQJ method. Here, 
|\&i(t)) and |^o(*))> which are defined in the total sys- 
tem Hilbert space Hs^Hs, are the possible values of the 
stochastic wave function of the combined GAW process. 
If the system part belonging to Us of the GAW stochastic 
wave function is in the same projective ray as the values 
of the stochastic wave function of the NMQJ method, we 
can conclude that both methods generate similar realiza- 
tions for the reduced system. We will also see that the 
deterministic evolutions of the stochastic wave functions 
for both processes are identical. 

This means that the PDPs of the two methods are the 
same in the following sense. The state space consists 
of the same set of projective rays in Us, stochastic wave 
functions evolve similarly between random jumps in both 



processes, and random jumps in both processes take place 
between the same two projective rays in Hs with equal 
rates. It is sufficient that the states belong to the same 
projective ray in Us since we are interested only in the 
dynamics of the reduced system. Moreover, the states 
in the same projective ray give equal contribution to the 
density matrix of the system since the complex phase of 
the state is not an observable. 

The summing of the GAW probability currents means 
that we lose the information to which mode the excitation 
from the system goes as the system decays. It is intuitive 
that the sum of the probability currents corresponds to 
the decay rate since the decay rate describes the total 
effect of the environment onto the system. However, it is 
important to note that for a given sign of the decay rate, 
there typically occurs probability flow components of the 
GAW realizations to both directions. 

In examples below, we will set the frequency dependent 
couplings to be real- valued and equal to = yjdvpk(vk), 

1 *Y A 2 

where Pkiy) = ^ ^ v _^ ')2+a 2 1S the spectral density, dv is 
the mode spacing, A is the spectral width, uj c is the po- 
sition of the peak in frequency space, and 70 defines the 
height of the peak. These parameters are also related to 
the time scales involved in the dynamics. We have 75 ~ 
7(7 1 , which is the time scale of the reduced system evolu- 
tion, and T£ ~ A -1 is the time scale of the environmental 
correlation functions. We can also compare our discrete 
TV-mode cases to the exact and numerical solutions ob- 
tained in the continuum limit J2k \dk\ 2 ~ * J dv pk{v). 

In the following, we make a detailed study for a two- 
level system (TLA) and a three level atom in a V- 
configuration ( V-system) . 



A. Two- level atom 

The Hamiltonian in the interaction picture is now 

N 

Hi = l Y / 9k(\g)(e\ale^ t - |e) (g\a k e- in "% (25) 
fc=i 

where Q k = v k ~ u eg . The state of the total system and 
the initial conditions are 

JV 

l*(*)) =(c 9 (t)\9) + c e (t)\e))\0)+J2ck(t)\g)\l k ), 



k=l 



Cfc(O) =0, 



(26) 



so that initially the modes of the environment are in a 
vacuum state. The Schrodinger equation and the inter- 
action picture Hamiltonian lead to the following system 
of first order differential equations for the amplitudes: 



,(*)=0, 



N 



fc=i 



c fc (t) = g k e tnkt c e (t). 



(27) 
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Probabilities to find zero or one photon in the environ- 



ment are from Eqs. (20) and &21 



p(o,t) = (nmv\*(t)) = \c g (t)\ 2 + \ce(t)\\ 



N 



p(M)H*(*)|ni|*(«)> = £Mt)| 5 



fe=i 



p(o,t). 

(28) 



In the GAW method, the combined property states, 
which are the two possible states that the stochastic wave 
function can take, are by using Eq. ( 22 ) 



l*o(*)) 



l*i(*)> 



c g (t)\g) + c e (t)\e) 



1 



1 x 



N 

^c k \g)\l k 



[0> = |0 o (t)>|O>, (29) 



^2c k \g)\l k ). 



(30) 



Here in the upper equation we use \<t>o(t)} to denote the 
reduced system part of the corresponding total system 
property state. From Eqs. (10) and (23) we get the com- 
bined probability current 



-2Re 



Cejt) 
Ce(t) 



\Ce{tt 



(31) 



The probabilities P(l,t) and P(0,t) satisfy Eq. (24) 



which can be easily calculated by using the Hamiltonian 
and the total state of the system and the environment, 
or the definitions P(l, t), P(0, t) and J\fi[t) [see the text 
below Eq.p4])]. We can define the transition rates by 
using Eqs. ||J and @. When Ji, (t) > 0, 

T ... . m , , \Ce(t)\ 2 



-2Re 

75, l (*) = 0, 
and when Ji${t) < 0, 
7i,o(i) = 0, 

Ce(t) 



c e (t) 

c e (t)} \c g (t)\ 2 + \c e (tW 



To ,i(t) = 2Rc 



\CeW 



C e (t)j 1 - |c ff (*)|= - |Ce(*)| ! 



(32) 



(33) 



Finally, the reduced density matrix can be obtained by 
taking the trace over the environment 

Ps (t) = its {W (t) | #0 (t)) {* o (t) I + W! (t) | *! (t)) <* 1 (t) | } 

= «;o(i)|^o(i))(^)(*)|+«'i(*)|ff)(fl|. 

Next we will study the TLA with the NMQJ method 
keeping in mind the previously derived results with the 
GAW method. The master equation for the TLA unrav- 
eled with the NMQJ method is 



d^ W 



A(t)(a.p s (t)a + - ±{ps(t),(r + v-}\ ,:>i4i 



where the decay rate A(t) and Lamb shift rate S(t) are [1 

Ce(t)' 



A(t) 
S(t) 



2Re 



21m 



c e (t) 

C e {t) 
Ce(t) 



(35) 



and the non-hermitian Hamiltonian giving the determin- 
istic evolution of the stochastic wave function is pQ 



Hes(t) = \ [S(t) ~ iA(t)} a+a. 



(36) 



All the amplitudes Cj(i) in Eq. (26) are solutions of the 



Schrodinger equation for the system and the environment 
with the Hamiltonian from Eq. (25). These amplitudes 



have the following connection to the normalized state vec- 
tors of the effective ensemble of the NMQJ method 

I , (t)) = c g {t)\g) + c e (t)\e) 
^\c e {t)\> + \c g {t)f 

\Mt)) = Iff), 
b g (0) =c g (0), b e (0) =ce(0). (37) 

Comparing these with the property state |\E , o(^)) of the 
GAW method in Eq. |29l, we can see that 



tr £ {|¥ (t)X*o(*)l} = \Mt))(Mt)\ = \Mt))(Mt)\, 

tre{mt)){*iV)\} = \9}{9\ = hM*)XlM*)l- (38) 

This shows that the reduced system part of the GAW 
realizations and the NMQJ realizations are identical. We 
are left with showing in detail that also the transition 
rates are the same. 

The reduced density matrix in NMQJ is 

ps(t) = P(\Mt))>t)\*a(t))(Mt)\ 

+ p(\Mt)),t)\Mt))(Mt)\- (39) 

When A(i) > we have transitions from \ip (t)) — s- 
\ipi(t)) and from Eq. (15) we obtain 



R lfi {t) = A(t)- 



|Ce(*)| 2 



C 9 {t)\ 2 + \Ce{tW 



(40) 



This is identical to 7i,o(i), when J\fi(t) > 0, see Eqs. (32) 
and (35). 



When A(f) < we have transitions from IV'i(i)} 
|^o(t)) and 



RoAt) = -T^yWh ice(<)|2 



,(i)| 2 + |c e (*)| S 



(41) 



Since ps(t) must be a positive operator we know that 
the decay rate A(t), and therefore also J\${t), must ini- 
tially be positive. Let us call t\ the time when A(t) 
turns negative for the first time. Now ps(t), when 
t < t\, generated by GAW and NMQJ must be the 
same, since from Eq. (38) we see that the states in 



(> 



the decomposition of ps(t) belong to the same projec- 
tive ray and R\ o (t) = 7i o (0 f° r i < ii ■ Therefore 
we have P(0,t) = \c g (t)\ 2 + \c e (t)\ 2 = P(\i> (t)),t) and 
P(l,i) = P(\ipi(t)),t). Now we can rewrite i?o,i(^) as 



Ro,i(t) = - 



P(0,t) \c e (t)\ 2 
P(l,t) {> P(0,t) 

\Ce(t)\ 2 



A(t) 



P(l,t)> 



(42) 



which is the same as %,i(t), when Jx,o(t) < 0, see 
Eqs. (33) and (35 1. It is also clear now that at t = t\ 



both J\ o(t) anchA(i) turn negative. 

Thus we have shown that we can derive the NMQJ re- 
sults from the GAW method for this system. This means 
(i) that we can obtain the decay rate in the master equa- 
tion ( 34 ) from the probability currents between the total 



system property states of the GAW method, and (ii) that 
the random state vector in %$ m both methods obtains 
its possible values from the same set of states, namely \g) 
and (we neglect the global phase since it plays no 
role here). 

In the first example we have chosen the parameters 
as \%{t)) = | e> 1 0> , time scale [t] = 1/A, 8 = 3A and 
7o = 0.8A. We use 180 environmental modes and a sta- 
tistical ensemble with 10 4 members. With the parame- 
ters mentioned above, the decay rate in the master equa- 



tion ( 34 1 is time dependent but always positive, thus cor- 



responding to the time-dependent Markovian case [15] , 
This also means that there are no reverse jumps in the 
NMQJ method in this parameter regime. However, as 
Fig. [I] shows, there are negative probability currents in 
the GAW method for specific modes or individual prop- 
erty states, while the total probability current between 
the system and the environment, J7i.o(£)j remains pos- 
itive indicating net current from the system to the en- 
vironment. This means that while there are individ- 
ual transitions from one photon to zero photon states in 
GAW, the number of transitions from zero photon states 
to one photon states is larger keeping the total probabil- 
ity current positive, which then matches the probability 
current obtained from NMQJ. 

In the second example we have chosen the parameters 
as in the first example except for 70 = 4A and 8 = —AX. 
The system is now in the non-Markovian regime display- 
ing also negative values for the decay rate. Figure[2]shows 
the individual probability currents for this case. The re- 
sults show that the region Vf. — uj c w (or « oj c ) 
gives the dominant contribution to the total probability 
current and it has also dominant negative contribution. 
As a consequence, the total current has negative periods, 
which is reflected in the negative regions for the decay 
rate, and thus the system is driven to the non-Markovian 
regime. 

In Fig. [3] we have plotted the decay rate which is cal- 
culated from the probability current components. Wc 
compare it to the exact decay rate calculated in the con- 
tinuum limit and see that the agreement of the curves 



is good. In the same figure we have also plotted the ex- 
act solution for the density matrix and compare it to the 
simulated ones, and we can see that the agreement of the 
curves is excellent. 




0.02 



0.015 



0.01 



0.005 



-0.005 



-0.01 



FIG. 1. Probability currents Ji k ,o(t) for TLA in the Marko- 
vian case. Initial state is |e}|0), [t] — 1/A, 8 — 3A, 70 = 0.8A 
and we use 180 environmental modes. When i ~ 1 we see 
that there occurs negative probability currents. 




FIG. 2. Probability currents Ji k ,o(t) for TLA. The parame- 
ters are as in Fig. [I] except that 70 = 4A and 8 = — 4A. We 
can identify the modes for which « uj c responsible for the 
non-Markovian effects. See the text for details. 



B. V-system 

The Hamiltonian for the V-system in the interaction 
picture is 



TV 



Hi =i^2,g k {\c)(a\a\e 



\c){b\a\e 



h.c), 



k=l 



(43) 
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FIG. 3. Decay rates and excited state population for TLA. 
The parameters are as in Fig. [2] 



where n k ,i = v k — ^i, i — o,,b, and we have denoted with 
| a) and |6) the two upper states and with |c) the ground 
state. Differential equations for the amplitudes obtained 
from Schrodinger equation are 



c b (t) 
c a (t) 



0. 



Cfc(t), 



JV 

fe=i 

JV 

fc=l 



(44) 



As we have seen in Sec. |IV A| the sum of the probability 
currents over all modes is related to the decay rate. By 



using the same procedure as in Sec. IV A it is possible 
to derive the following equations for the probabilities 



At 
_d 

d7 



P(0,t) = -Ji,o(*), 
P(l,t) = J 1>0 (t). 



(45) 



The reduced system dynamics corresponding to these 
equations is given by a non-secular master equation 
which is not in general compatible with the form given 
in Eq. (12) used as a starting point for the NMQJ. 



To find the connection between GAW and NMQJ in 
this system, we approximate the exact non-secular dy- 
namics by de-coupling the evolution of the coherences 
and populations |28j . Eventually this means that the 
emission of the photon can be associated to one of the 
two decay channels and we write the Hamiltonian as 



H, 



N 

fe=i 



g k (\c){a\a{e 



1 i 1 1, 



\c)(b\ble 



+ h.c), 



(46) 



where we have introduced new environmental modes de- 
scribed by operators b k . This means that we can identify 
from which decay channel the photon originated, which 
prevents the occurrence of quantum beats [35] . 

The differential equations for the amplitudes are then 



c c (t) = 0, 
c b (t) 

C a (t) 



JV 

E 

fc=i 

JV 



9ke 



■Eff* <r<n * ,B * c £(*)' 

k=l 

c b k (t) = g k c b (t)e m ^\ 



(47) 



and these equations are a good approximation for 
Eq. (44) in certain parameter regions. 



We assume that initially the environmental modes are 
empty. Then we can give the total state of the system 
and the environment as 

\V(t)) =c c (t)\c)\0) a \0) b + c a (t)\0) a \0) b + c b (t)\0) a \0) b 



JV 



^c%(t)\c)\l k ) a \0) b + ^cUt)\c)\0) a \l k ) 



k=l 



k=l 



(48) 



From now on we drop the subscripts referring to different 
Hilbert spaces. We want to know the probabilities to find 
a photon in the environment and we want to identify 
which of the excited states has decayed. Therefore it is 
natural to use the following operators 

n = / s ®|o)(o|®|o)(o|, 

JV 



ni,o = 5^/5®|ifc)<ifc|c 

k=l 

JV 

ni, 6 = ^7s®|o)(o|®|i fe )(ifc 



(49) 



k=l 



The probability to have one photon in the environment 
which has been created when the excited state i decayed, 
is P l (l,t), where i = a, b, and the probability to have 
zero photons in the environment is P(0,t). Following a 
similar procedure as for the earlier presented TLA case, 
we can calculate them as 



JV 



P i (i,t) = <*(t)|n 1)< |*(t)) = x;i4(*)l 2 . 



k = l 



p(o,t) = <*(t)|n |*(t)) = | Cc (t)| 2 + \c b {t)\ 2 + | Ca (i)| 2 . 

(50) 



Subsequently, the combined property states are 

c c (i)|c)|0)|0)+c o (t)|a)|0)|0) + c b (t)|6)|0)|0) 



l*o(t)) 

l*l,a(*)> 
l*l,6(*)) 



(51) 



Ef=i^W|c)|i fe )|o) 



ZLi4(t)\c)\o)\W 



The differential equations for probabilities P (0,t ) and 
P l (l, t) can be calculated with the help of Eqs. Q, @ 
and d48l. We obtain 



dt 



P(0,t) = -Mt) - jUt), 



dt 



dt 



P a (l,t) = J^t), 

p b (i,t)=jut), 



(52) 



where the combined probability currents J\o{t) and 
J[ o(t) tell how much probability is flowing from the sys- 
tem to the environment in each channel separately. The 
combined probability currents are now 

J^ (i) = -2Re|^||||c l (0| 2 , (53) 

where i = a,b. We can define transition rates as in 
Eqs. ^ and ^ separately for each decay path since we 
can partition the combined probability current into two 
independent parts. They are, when J{a{t) > 0, 



rut) 



P(0,t) 



-2Re 



Ci(t)j p(o,ty 



7?,i(*) = 0, 
and when Jl {t) < 0, 

V. (t) = o, 



(54) 



VAt) 



>(l,t) 



2Re 



Ci(t) \ \Ci(t)\< 



a{t) J P(i,ty 



(55) 



namics under secular approximation is 

^5(«) = -i[^ (t)|a)(a| )Ps (t)]-i[^ 6 (t)|6)<6|, /9s (t)] 
+ A a (t)(|c><a|p s (t)|a)< C |-^{ Ps (t),|a><a|} 



A b (t)(\c)(b\p s (t)\b)(c\--{ Ps (t),\b)(b\} 



(57) 



The total state of the system and the environment has 
been given in Eq. ( 48 1 , and by tracing out the environ- 



mental degrees of freedom and taking the time derivative 
of the expression we can identify the Lamb 

shifts and the decay rates to be 



A<(t) 
Si(t) 



2Re 
-21m 



MO 

Ci(*) 
Cj(*) 
Ci(t) 



(58) 



where i = a, 6. The non-Hermitian Hamiltonian for 
NMQJ in this system is 



#cffW=E^(i)-*A 4 (t))|z)< 



(59) 



where again i — a,b. We can give the deterministically 
evolving state of the NMQJ process and the initial con- 
dition as 

\Mt)) = d c (t)\c) + da(t)\a) + d b (t)\b), (60) 
d i (0)=c i (0), 

where j = c,a,b and Cj(t) are probability amplitudes 
from Eq. ( 48 1 . By solving the time evolution given by 



the Hamiltonian of Eq. (59), we see that c c (t) = d c (i), 



c a (t) — d a (t), and c b (t) = db(t). In NMQJ, the realiza- 
tions of the process are normalized and therefore we can 
write 

= C (t)\c) + c a (t)\a) + c b (t)\b) 
V|c c (t)| 2 + |c h W| 2 + k(t)| 2 
\Mt)) = \c) 

The reduced density matrix of the NMQJ process is then 



where i — a,b. 

The reduced density matrix generated by the GAW 
method is now 

Ps(t) =trg{«;o(t)|*o(t))<*o(t)|+t«l l a(t)|*l l a(*))<*l J a(t)| 
+ «»1,6(*)I*1,6W)(*1,6(*)|}- (56) 

Next we will study the NMQJ method for this system. 
The master equation describing the reduced system dy- 



Ps (t) = p(\Mt)),t)\Mt))(Mt)\ 
+ p(\Mt)),t)\Mt))(Mt)\- 



(62) 



As in Eqs. (15) and (18), and using Eqs. (50), we can 



write the transition rates, when Aj(t) > 0, 

h(tW 



R\o(t) = Mt) 
R} )1 (t) = 0, 



p(o,ty 



(63) 
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and when Aj(i) < 0, 
R\ o (t) = 0, 



P(\Mt)),t) 



Ai(t) 



\Ci(t)\ 



p(IVi (*)},*) w P(o,t) 



(64) 



where i = a,b. By using Eqs. (51 ) and (161]), we obtain the 



connection between the GAW and NMQJ state vectors 



tr £ {|*o(*)X*o(t)|} = |^(*)X^)(*)l, 
trg{|* lii (t))<* M (t)|} = ^i(t))<Vi(*)l, 



(65) 



where i — a,b. 

This means that the system Hilbert space part of the 
possible realizations of the combined GAW process and 
the NMQJ process with the same index belong to the 
the same projective ray in T-Lg. We also see that there 
is redundancy in tr^ {|\E f i..j(t)){ 1 I f i.i(t)|} since both states 
when i = a, 6 belong to the same projective ray in %s- 
It means that we can combine w\^{t) + w\ ta {t) — w\{t) 



in Eq. (56) and that both transition rates T a and T 



induce jumps between the same two projective rays but 
the rates of the jumps are generally different. 

We assume that from some initial time to to t\ the rates 



A a (io) an d A;,(i) are positive. Then from Eqs. (54), (58) 
and pj) we see that R\ (t) = T-} {t) and i? 01 (t) = Tjjt). 



This implies that the density matrices generated by GAW 
and NMQJ are the same at least to t\ when at least one 
of the decay rates or collective probability currents turns 
negative. Since NMQJ and GAW both have normal- 
ized realizations we can deduce that WQ(t) — P(0,t) — 
P(\Mt)),t) and Wl (t) = P(l,t) = P(\Mt)),t) when 

t€[t ,h]. 

Now when negative currents (J~l (i) < 0) or decay 
rates (Aj(t) < 0) emerge when t > t\, the transition 
rates are 7q i(t) — R l and Ti (t) = R\ (t) since we 
have P(0,tj = P(\Mt)),t) and'P(l,*) =' P(W>i(t)), t) 
(the calculation is the same as we did in Eq.(|42j>). Thus 
we have shown that the GAW process for the combined 
property state is an equivalent process to NMQJ in Us 
in the sense we defined in the beginning of Sec. |IV| 

Next we study a numerical example where the ini- 
tial state is written in non-secular case as ^(i)) = 
^(|a)s|0)£ + |6)s|0)£) and under secular approxima- 
tion, where each channel has independent environment, 
as !*(*)) = ^(|a> s |0) £ JO) £t + \b) s \0) £a \0) £b ). We are 
written here the different Hilbert spaces explicitly for 
clarity but from now on we omit this for compactness 
of notation. Other parameters are defined as [t] = 1/A, 
7o = 4A, 5 a = 3A, 5b = — 3A and we have used 240 envi- 
ronmental modes. 

We start with the non-secular case. Figure |4] shows the 
corresponding probability currents between the property 
states of the GAW method. The interference of prob- 
ability currents is visible. In the secular case Fig. [5] 
shows probability currents under the secular approxima- 
tion which decouples the two excited states. Probabil- 
ity currents can not interfere because each excited state 



interacts with its separate environment. The compari- 
son between the density matrices for the two cases are 
shown in Fig. [6j In the same figure we can also see that 
under secular approximation the reduced dynamics are 
governed by master equation (57). We can clearly see 



the effect of the interference of the probability currents 
to the reduced system dynamics in the non-secular case. 

In Fig. [7] we show the effect of the secular approxima- 
tion to the combined probability current from the sys- 
tem to the environment, ie. the difference of »7i,o(t) an( i 
^oM+^oW- There are fast oscillations in J\fl(t) and 
it can even be negative when J^ {t) +J\ a{t) is positive. 
In Fig. [7] we compare the secular approximation decay 
rate calculated from GAW to the TCL2 [Tj decay rate 
and we see that the match is very good. 



/ / / / 1 j 1 1 / j / \ / 




1 



0.05 
0.04 
0.03 
0.02 
0.01 


-0.01 
-0.02 
-0.03 
-0.04 
-0.05 



FIG. 4. Probability currents in the non-secular case. The 
initial state is 4g (|o)|0) + 16)|0)), [t] = 1/A, 70 = 4A, 5 a = 3A, 
5b = — 3A, and we have used 240 environmental modes. 




0.025 




FIG. 5. Probability currents in the secular case. The pa- 
rameters are as in Fig. HI but the initial state is ^j(|i}|0}|0} + 
|6)|0)|0». 
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Time Time 



FIG. 6. Density matrix elements for the V-system The pa- 
rameters are as in Figs. [4] and [5] For the chosen parameter 
values p aa (t) = Pbb(t). 



FIG. 7. Top: Sum of probability currents from the system 
to the environment. Solid red line is the non-secular case 
where we can not distinguish different decay channels. Blue 
circles are the secular case where we have two different decay 
channels. Bottom: Decay rates calculated from GAW method 
for the secular case (solid blue line) and decay rate of TCL2 
master equation (red circles). Parameters are as in Figs. [4] 
and[U 



V. CONCLUSIONS 



We have studied the non-Markovian dynamics of sim- 
ple quantum optical systems by means of two jumplike 
unravelings. The GAW method uses piecewise determin- 
istic realizations within the Hilbert space of the total 
system while the NMQJ method exploits piecewise de- 
terministic realizations within the Hilbert space of the 
reduced system. Our analysis shows that there exists 
a connection between the two methods. In particular, 
we have demonstrated that the reduced system part of 
the property states of the GAW are identical for the 
NMQJ state vectors in the considered cases. Moreover, 
the summation over the probability currents appearing in 
the GAW formalism are directly connected to the decay 
rates of the time-local master equations and hence to the 
rates of jumps in the NMQJ method. While there exists 
quite a large variety of Monte Carlo methods for non- 
Markovian systems [TlTH3l I27H37] . both jump and dif- 
fusion type, generally the connections between the meth- 
ods have not yet been extensively investigated apart of a 
few studies J35j [44] . We expect that the results presented 
here stimulate further research in this area leading to im- 
proved insight to the often complex quantum dynamics of 
non-Markovian systems. Moreover, analyzing the prob- 
ability currents in similar manner as treated here, may 
lead to further understanding of the information flow be- 
tween the system and the environment, a topic which is 
currently vividly discussed in the context of open quan- 
tum systems. 
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